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ABSTRACT 

The  stability  of  incompressible,  gravity  free,  immiscible  two  phase  flow 
to  viscous  fingering  is  studied  numerically.  In  contrast  to  numerical  studies  in 
the  literature,  but  in  good  qualitative  agreement  with  laboratory  studies,  we 
find  a  large  parameter  regime  for  which  the  displacement  is  strongly  unstable. 
The  pjirameter  space  examined  includes  mobility  ratio,  capillary  pressure 
strength,  size  of  perttirbation,  and  heterogeneity  of  the  medium. 
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1.  IntrodDction 

Two  phase  miscible  and  immiscible  displacement  processes  occur  in  reservoir  operations 
and  in  laboratory  core  floods.  The  stability  of  miscible  and  immiscible  flood  has  been  stu- 
died by  phenomenological  models, l-^. 3.^  numerical  simulators.^-^  and  by  linearized  stability 
analyses. I'^'S'^  In  general,  phenomenological  models  assume  simplified  flow  descriptions  so 
that  other  effects  may  be  studied  in  more  detail.  Numerical  simulation  is  often  used  but  usu- 
ally at  a  cost  of  large  numerical  diffusion  effects.  Linearized  stability  analyses  provide  infor- 
mation on  the  inception  of  instability  for  small  perturbations  but  give  no  information  on  the 
effects  of  perturbation  size  or  the  eventual  form  of  fingered  profiles. 

The  purpose  of  this  study  is  to  perform  a  fundamental  examination  of  the  stability  of 
immiscible  flood.  The  usual  relative  permeability  and  capillary  pressure  functions  are  used  to 
describe  the  local  flow;  no  effective  (average)  "shape  factors''^-^  are  assumed.  Care  is  taken 
to  ensure  that  the  physical  effects  (capillarity,  heterogeneity)  dominate  the  numerical  disper- 
sion present  in  the  calculation. 

The  instability  of  immiscible  flow  is  studied  instead  of  miscible  flow  for  two  reasons. 
First,  it  is  more  difficult  to  model  numerically  than  miscible  flood;  the  Buckley-Leverett  mix- 
ing zone  leads  to  enhanced  stability  compared  to  the  corresponding  miscible  flood.  Second, 
in  addition  to  viscous  and  dispersive  forces,  miscible  flood  behavior  is  further  complicated  by 
the  phase  behavior  of  components,  non-trivial  chemical  interactions  between  components  and 
rock,  and  even  regions  of  immiscible  flow  (CO2  flood).  In  immiscible  processes,  the  instabil- 
ity is  dominated  by  flow  physics  (interplay  of  viscous  and  capillary  effects) ;  in  many  instances 
the  effects  of  chemistry  are  of  secondary  importance.  Therefore,  the  model  (1)  is  expected 
to  describe  a  realistic  flood  without  extensive  simplification  of  the  process. 

Laboratory  studies  of  flow  through  core  samples  exhibit  a  variety  of  fingering 
phenomena. i-'  Linearized  analyses  of  two  phase,  incompressible,  immiscible  flow^-^  reveal 
regimes  of  definite  instability  for  the  growth  of  viscous  fingers.  Published  numerical  studies^ 
in  the  non  Unear  range  show  little  evidence  of  instability  for  even  moderately  strong  adverse 
mobility  ratios.  The  results  of  our  study  show  that  understanding  of  the  interplay  of  length 
scales  present  in  the  numerical  study  is  essential  for  a  correct  numerical  calculation  of  the 
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stability  of  viscous  fingering . 

The  problem  can  be  chararterized  in  terms  of  four  dimensionless  length  scales: 

1)  ^cap  ™  ^cap  I  ^  characteristic  of  the  capillary  strength. 

2)  ^htx  ~  ^hei  I  ^  characteristic  of  either  the  heterogeneity  in  the  medium  (or  characteristic 
of  the  length  scale  of  perturbation(s)  inserted  in  the  initial  data  if  the  medium  is  homo- 
geneous). 

3)  ^  ■  L  /  L  •■  max(Ajr,  \y)  /  L  where  Ajc  and  tiy  are  the  mesh  spadngs  (assumed  reg- 
ular) of  the  grid. 

^)  ^disp  ™  ^diip  ^  ^  characteristic  of  a  minimum  'fiducial  wavelength*  for  the  numerical 
scheme  due  to  inaccuracies  arising  from  tnmcation  error  of  the  finite  difference  scheme 

employed. 

In  the  above  L  is  a  reference  length  scale  which  we  shall  assume  to  be  Z,  =  min(Lj,  L^) 
where  L,  and  L^.  are  the  sizes  of  the  computational  domain. 

The  dispersion  present  in  the  numerical  scheme  is  present  in  two  forms,  a  purely  disper- 
sive piece,  and  a  diffusive  part,  corresponding  respectively  to  the  real  and  imaginary  parts  of 
the  dispersion  relation  ^(Jl:)  (w  "  frequency,  k  "  wavenumber)  for  the  difference  scheme. 
The  dispersion  relation  can  easily  be  derived  for  a  numerical  scheme  (such  as  used  here) 
when  it  approximates  a  linear  partial  differential  equation  (PDE)  such  as 

u,  +  a  u,  =  b  u„ 

where  a  and  b  are  constants,  with  the  usual  results  that  the  dispersion  relation  (both  real  and 
imaginary  parts)  for  the  numerical  method  matches  the  dispersion  relation  for  the  PDE  to 
within  reasonable  accuracy  for  wavelengths  longer  than  some  minimum  (fiducial)  wavelength 
Xyj^.  The  value  of  X^^  is  dependent  on  the  ratio  a/b.  For  non-linear  PDE's  such  as  investi- 
gated here,  the  dispersion  relation  is  not  amenable  to  calculation  and  the  results  from  linear 
analyses  have  to  be  assumed  to  hold  locally  and  used  as  guidelines.  For  the  second  order 
accurate  method  employed  here,  investigations  on  the  above  PDE  show  such  that  X^  ^  5  L^ 
to  be  an  expected  reasonable  fiducial  length  for  the  parameter  range  discussed  in  this  study. 
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The  problem  investigated  here  has,  effectively,  a  smallest  physical  length  scale  deter- 
mined by 

Care  was  taken  to  ensure  that  the  relation 

held  for  the  calculations  done  here/  However,  though  due  to  computational  memory/time 
limits,  there  are  local  areas  in  some  of  the  calculations  where  R^^^  &  3/?^.  Some,  isolated, 
local  distortion  of  the  calculation  results  (c.f.  Fig.  8c),  but  does  not  affect  the  main  conclu- 
sions of  the  calculation. 

In  addition  to  the  above  restrictions  however,  we  find  that  if  the  length  scale  over  which 
the  heterogeneity  varies,  /?;,,,,  is  less  than  the  effective  diffusion  length  scale, 
max(R^,jp,  R^^p),  the  stability  properties  of  the  flow  are  the  same  as  if  the  medium  was 
homogeneous.  (See  §  4.2  .)  This  is  a  physical  consequence  of  the  diffusion  term  and  reflects 
the  enhanced  stability  emparted  to  the  flow.  This  we  propose  is  the  main  reason  for  failure 
of  past  efforts^  to  resolve  the  instability  of  immiscible  flow  numerically. 

2.  The  Displacement  Process 

We  consider  the  immiscible  displacement  determined  by  the  following  system, 

*  IT  "^  VV,  =  0  ,  (la) 

V,  =  -his J  VP,  ,         («•  -  o{il),  w{ater))  ;  (lb) 

where  S^  +  S^  =  1,  and  \,  =  K(x,  y)  k^,(Si)  I  ji,  is  the  mobility  for  component  i.  The  rock 
permeability  K{x,  y)  is  taken  to  be  a  scalar  rather  than  a  tensor  field.  Defining  a  mobility 
weighted  average  pressure 


'  The  so-called  'grid  orientation'  problem  is  a  well  known  example  of  the  effects  of  numerical  trunca- 
tion error  present  in  full  scale  reservoir  simulators  where  the  length  scale  of  a  grid  block  is  almost  al- 
ways much  larger  than  the  smallest  physical  length  scale  associated  with  the  problem.  Indeed,  for 
Buckley-Leverert  or  plug  flows  where  the  oil-water  interface  is  effectively  of  zero  width,  numerical 
schemes  must  be  employed  w-hich  either  have  the  same  local  dispersion  charaaeristics  as  the  partial  dif- 
ferential eauation  being  solved'-"  or  that  have  a  dispersion  charaaeristic  which  has  controlled  angular 
behaviour." 
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^  _   (X,  P^  +  K  Pq)  ^  (2) 

A 

where  \  —  X^  +  X^,  and  the  capillary  pressure 

Pc^'Po-  p.,  (3) 

system  (1)  can  be  written 

«}>  -^  +  V  •  (/(5J  V  )  =  V  •  (  g{SJ  V  S,  )  ,  (4a) 

V  •  V  =  0,  (4b) 


where 


V  =  -X(5J  {   V  P  +  P,(;e,  5J  /'(5J  V  5,   }  ,  (5) 


,^  ,_    X.,(5J  X,(5J 


f     dPXSJ 


ds 


(7) 


In  (4),  S^  and  /*  are  taken  as  the  fundamental  variables  and  P^  is  assumed  to  be  a  known 
function  of  S„  and  j?.  The  position  dependence  occurs  in  heterogeneous  media.  The  mobility 
weighted  average  pressure  (2)  has  the  virtue  over  the  unweighted  average  (P„  +  P^)  /  2  of 
being  finite  for  all  S„. 

K  P,(S^)  ^  0,  P„  -  P„,  and  (4)  reduces  to  the  familiar  equations  of  Buckley-Leverett 
flow.  However,  (4)  with  P^  "  0,  is  not  a  well-posed  problem  in  the  unstable  flow  regime; 
small  wavelength  disturbances  grow  at  a  rate  inversely  proportional  to  their  wavelength.  One 
must  require  physical  dispersion  (eg.  capillarity)  dominate  over  numerical  dispersion  to 
abtain  a  physically  meaningful,  well-posed  problem. 

3.   The  Namerical  Approximation 

Equation  (4a)  is  degenerate  parabolic;  the  effects  of  changes  in  saturation  in  one  area 
propagate  to  surrounding  areas  with  a  finite  speed.  Equation  (4b)  is  elliptic  in  the  pressure 
P;   all. areas  affect  each  other  simultaneously  -  the  propagation  speed  is  infinite.    Because  of 
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the  difference  in  time  scales  between  these  two  equations,  they  can,  with  a  great  deal  of  con- 
fidence, be  solved  sequentially.  One  cycle  consisting  of  solution  of  (4b)  followed  by  (4a) 
shall  be  referred  to  as  one  iteration  (or  timestep).  The  main  criterion  for  choice  of  the 
numerical  scheme  used  to  approximate  (4a)  and  (4b)  is  that  it  introduce  no  numerical  diffu- 
sion to  mask  the  physical  diffusion  produced  by  the  capillarity.  This  is  achieved  by  approxi- 
mating all  derivatives  by  central  differences. 

With  J,  j  and  n  denoting  respectively  the  x,  y  and  t  indices,  the  numerical  scheme  used 
for  (4a)  consists  of  advection  terms  discretized  by  a  leap-frog  stencil,  the  viscosity  term  by  a 
Crank-Nicholson  stencil: 

2  Af  Ax  A:y 

=  — ^ — Ip"      a  iu"*^  ■¥  uV^)  +  J?"      A  (u""^^  -t-  u"~^n 
4  Ax^  1  '  •■'     '    '  ••'         '  •■/  '''i  "^'^"rj       "r,]>j 

- — I?"      A  (u"^^  +  i/""^)  -^  b"      a  (m""^^  -t-  m''~^)1  , 


4Ay^ 


(8) 


where 


8,*  J  ■  8 


81- J  -  8 


2 


"f.y  ^  "I- 1, J 


A,  u'l^j  -  uf^^j  -  ulj  , 

Equation  (4b)  has  the  form 

V  •  X(5)  VP  =  V    his)  V5  , 
and  was  approximated  by  a  Crank-Nicholson  stencil: 


(9) 


(10) 
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The  difference  scheme  is  second  order  accurate  in  space  and  time.  The  linear  system  of 
equations  generated  on  the  numerical  grid  for  each  equation  at  each  timestep  were  inverted 
using  a  multigrid  algorithm^  of  relaxation  and  injection  through  a  sequence  of  coarser  grids 
which  is  based  on  a  variational  formulation  of  the  equations  (4). 

The  calculations  were  performed  in  a  square,  length  L,  width  L.  The  boundary  condi- 
tions imposed  were  (5  is  defined  in  (11)) 

5=1,    at^^O,  5  =  0,    aXy  =  L  , 

—  =  0        at  X  =  0,  Z,  , 

dx 

P  =  Po,    at  >-  =  0,  P  =  0,    aty  =  L  ,  (11) 

—  =  0        &tx  =  0,L  , 

dx 

uix,y,  f  =  0)  =  UQix,y)  . 

The  calculation  neglects  the  initial  displacement  flow  into  the  block  during  which  the  water 
saturation  at  >  =  0  rises  from  5^^  to  1  -  S^^.  It  is  not  expected  that  the  initial  flow  into  the 
block  will  affect  the  stability  properties  of  the  displacement  and  hence  this  initial  transient 
flow  is  neglected.   Further,  this  choice  simplifies  the  saturation  boundary  condition. 

The  constant  pressure  drop  boundary  condition  contrasts  with  the  (perhaps  more  com- 
mon) experimental  condition  of  maintaining  constant  flow  rate.  Constant  flow  rate  calcula- 
tions can  be  included  but  are  more  complicated  due  to  the  non-uniform  flow  at  y  =  0. 
Again,  the  stability  properties  of  the  system  are  not  expected  to  be  qualitatively  different  for 
the  two  types  of  boundary  conditions. 

4.   Analysis  of  Case  Stadies 

In  the  absence  of  capillarity,  the  Buckley-Leverett  mixing  zone  terminates  in  a  sharp 
oil-water  interface,  with  a  unique  water  saturation  5*  just  behind  the  front.  If  5^^  denotes 
the  connate  water  saturation  ahead  of  the  front,  the  frontal  mobility  ratio, 
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is  a  measure  of  the  potential  instability  of  the  front  to  viscous  fingering.  The  'far  field' 
mobility  ratio, 

M  -   '^'"V  ,  (13) 

is  a  measure  of  the  potential  instability  of  the  front  only  in  the  case  of  two  regions  of  one- 
phase  flow,  i.e.,  plug  flow.  The  importance  of  M  and  M'  is  verified  by  linearized  stability 
analyses  of  plug  flow  and  of  Buckley-Leverett  flow.  Although  genuine  two-phase  flow  is 
usually  encountered,  M  is  commonly  used  instead  of  A/*  as  an  index  of  instability  because  it 
is  more  easily  measured  in  routine  laboratory  work.  In  general,  M  >  M',  and  thus  there  is  a 
region  of  parameter  space  for  which  M'  indicates  stability  while  M  indicates  instability. 

When  capillary  forces  are  included,  the  sharp  oil-water  front  disappears,  replaced  by  a 
region  of  smooth  rarefaction  from  5*  to  S^  over  a  length  scale  L^p  characteristic  of  the  local 
capillary  strength.  The  mixing  zone  is  thus  increased  over  the  Buckley-Leverett  case,  and  the 
system  is  expected  to  have  further  enhanced  stability  against  viscous  fingering.  As  L^^ 
increases,  both  M  and  M'  become  poorer  indicators  for  the  stability  of  the  displacement. 
Assimiing  quadratic  permeability  functions, 

k^  =  k^s\       k^  =  k,„{i-sy,  (14) 


where 


1   ~   ^onv   ~   ^wc 


(15) 


the  far-field  and  frontal  mobility  ratios  for  Buckley-Leverett  flow  are  given  by 


M' 


M  =  ^^,  (16) 


=  2(1  -  (1  -I- A/)-^')  .  (17) 


From  (16)  and  (17),  0  ^  M  <  «  whereas  0  ^  Af*  <  2,  reflecting  the  increased  stability  of 
immiscible  flow. 

It  is  expectedl3.i4  ^^^  ^  variety  of  factors  affect  the  stability  properties  of  two  phase, 
immiscible,  incompressible  displacement.  These  include 
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mobility  ratio  of  the  fluids  present 

strength  of  capillary  forces 

size  of  existing  perturbations  in  the  flow 

heterogeneity  of  the  reservoir 

injection  rate  of  displacing  fluid 

presence  of  injection  and  production  wells. 

The  first  four  factors  are  examined  in  this  study.   Because  of  our  boundary  conditions  the  last 
two  were  not  studied. 

As  it  is  computationally  impractical  to  pursue  a  full  mapping  of  the  parameter  space, 
our  calculations  present  projections  of  the  'stability  manifold'  along  a  portion  of  each  of  the 
pcirameter  axes  studied  and  show  the  development  towards  instability  for  each. 

All  computations  were  performed  in  a  square  (in  the  x,  y  plane)  slice  of  rock  of  poros- 
ity 0.18  and  average  rock  permeability  K  =  106  mD. 

The  pressure  drop  maintained  across  the  block  (in  the  y  direction)  is  such  that  if  the 
slice  had  a  cross-sectional  area  of  15.9  cm  (in  the  x,  z  plane),  the  flow  rate  would  be  0.00398 
cm^  /  sec.  at  f  =  0.  The  relative  permeability  functions  were  modelled  by  the  quadratic 
forms  (14),  and  the  capillary  function  represented  by 

P,  =  P^{x,y)]iiS    (atm).  (18) 

These  functions  and  constants  are  characteristic  of  sandstone.  ^^ 

Two  measurements  were  used  in  analysing  perturbation  (finger)  growth.  Figure  1  illus- 
trates the  measurement  of  the  height  (H)  find  full  width  at  half  maximum  (P^^HM)  for  a 
finger  present  in  a  particiilar  saturation  contour  5  =  5„  at  some  time  f  &  0  in  the  flow.  Sta- 
bility of  a  contour  level  is  examined  by  plotting  H  and  the  ratio  H  /  FWHM  for  a  finger  as  a 
function  of  time.   In  all  cases  only  the  contours  5  =  0.8,  0.6,  0.4,  0.2,  0.0  were  analyzed. 
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4.1.   Homogeneous  Rock  Studies 

For  homogeneous  rock,  {K(x,  y)  =  const,  ^0(0:,  >-)  =  const)  the  initial  data  Uq(x,  y) 
consisted  of  a  finger,  imposed  on  a  one-dimensional-flow  profile  lioCx)-  While  the  one 
dimensional  profile  could  have  been  obtained  from  one  dimensional  calculations  for  which 
the  water  saturation  at  ^  =  0  builds  up  from  S^^  to  1  —  Sg^w^  we  chose  instead  to  initialize 
MqCv)  wi^  3  quadratic  shape.  (See  Fig.  2c  .)  This  also  simplified  the  choice  of  Ijow  to  impose 
the  fingered  perturbation;  it  was  given  a  quadratic  variation  in  x,  y.  (See  Fig.  2a  .)  Conse- 
quently, our  calculation  includes  some  initial  transients  as  the  flow  adjusts.  We  emphasize 
again,  that  these  transients  are  not  expected  to  significantly  affect  the  stability  properties  of 
the  displacement  process;  to  whit:  if  a  certain  choice  of  parameters  in  our  calculation  indi- 
cates an  instability  of  a  certain  growth  rate,  the  instability  with  the  same  growth  rate  results 
whether  or  not  the  transients  are  included.  Figure  2  displays  a  selection  of  plots  showing  the 
perturbation  development  for  a  particular  nin  (run  3  of  §4.1.3).  Shown  is  the  saturation 
development  for  four  equally  spaced  timesteps  (Fig.  2a),  the  final  pressure,  saturation  and 
velocity  fields  in  the  run  (Fig.  2b),  and  the  cross  sectional  pressure  and  saturation  profiles 
along  the  mid-line  x  /  L  =  0.5  for  5  equally  spaced  timesteps  (Fig.  2c)  .  From  finger  growth 
measurements,  the  perturbation  shown  here  is  seen  to  be  slightly  unstable. 

4.1.1.  Stability  Dependence  on  MobUity  Ratio 

A  moderately  large  finger  was  initialized  on  a  one  dimensional  backgroimd  flow.  For 
the  S  =  0.8,  0.6,  0.4,  0.2,  0.0  contours,  the  size  of  the  perturbation  introduced  is  given 
below. 

sat  H  /  L  FWHM  /  L 

0.8  0.15  0.23 

0.6  0.21  0.23 

0.4  0.26  0.23 

0.2  0.30  0.23 

0.0  0.32  0.23 

Measurements  of  the  finger  development  on  these  contours  were  taken  at  constant  time  inter- 
vals during  the  run  for  three  different  mobility  ratios: 


-li- 
ning  M 

1  0.722       1.45 

2  1.0  3.0 

3  1.397      10.0 


In  all  cases  the  mobility  ratio  was  determined  by  choosing  m-w  =  1.  V^o  -  10,  kg„  =  1  and 
varying  Jk^  appropriately.  The  capillary  force  was  set  by  using  Pq  =  0.005  atm.  giving  a 
measured  capillary  length  scale  R^^^  ~  0.05  . 

The  results  of  finger  height  versus  injected  pore  volume  for  the  analyzed  contours  are 
presented  in  Fig.  3a.  As  expected,  the  flow  is  much  more  unstable  for  run  3  than  for  run  1  . 
All  contour  levels  analyzed  are  unstable  for  runs  2  and  3,  whereas  only  the  0.8  and  0.6  con- 
tours are  imstable  in  run  1."  We  note  that  the  growth  rate  of  the  height  (after  initial  tran- 
sients have  disappeared)  is,  to  a  first  approximation,  linear.  This  is  in  accord  with  results 
obtained  for  leirge  scale  perturbations  in  Buckley-Leverett  flow.  13 

Fig.  3b  plots  the  ratio  H  /  FWHM  versus  injected  pore  volume  for  the  same  three  nms. 

4.1.2.  Stability  Dependence  on  Initial  Finger  Size 

In  a  second  series  of  calculations,  the  mobility  ratio  and  capillary  pressure  constant  were 
held  constant  at  M'  =  1.397,  Pq  =  0.005  (R^p  ~  0.05)  and  the  size  of  the  perturbation  in  the 
initial  data  varied.   For  the  5  =  0  contour,  the  initial  perturbation  sizes  were^ 

run  H/Z.  FWHM/L 

1  0.09  0.07 

2  0.19  0.14 

3  0.32  0.23 

4  0.54  0.23 

The  results  are  plotted  in  Fig.  4a  and  b.    The  flow  becomes  more  stable  as  the  initial 


'  There  is  some  question  as  to  whether  the  initial  transient  flow  for  the  0.8  and  0.6  saturation  contours 
has  terminated  in  this  run  or  whether  it  is  stiL  adjusting  from  the  initial  conditions.  These  contours  are 
initially  isolated  from  the  frontal  flow  and  in  flows  of  borderline  stability  can  be  expected  to  remain 
t^us  isolated  for  a  longer  period  of  time. 

'  The  heights  given  in  column  two  of  the  table  are  numerical  approximations  to  0.1,  0.2,  0.33  and  0.55 
which  arise  due  to  the  finite  size  of  the  grid  spacing.  .  c  . 
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perturbation  size  decreases. 

4.1.3.  Stability  Dependence  on  Capillary  Pressure 

For  a  third  series  of  runs,  the  mobility  ratio  and  initial  perturbation  size  were  held  con- 
stant (A/*  =  1.397,  initial  perturbation  as  in  Fig.  3)  and  the  capillary  pressure  constant  P^, 
varied: 


run 

P, 

1 

0.005 

2 

0.01 

3 

0.02 

4 

0.05 

5 

0.1 

The  results  are  shown  in  Fig.  5.   Sufficient  increase  in  capillarity  is  seen  to  stabilize  the  flow. 

4.1.4.   Comparison  with  No  Capillarity 

The  front  tracking  code  of  Glimm,  McBryan,  et.  al.,10  was  used  to  provide  a  com- 
parison to  instability  growth  for  the  case  of  zero  capillarity  (Buckley-Leverett)  flow.  For 
these  calculations,  a  finger  was  initialized  in  a  homogeneous  medium  such  that  the  contour 
level  5*  at  r  =  0  occurred  at  the  same  position  as  it  would  in  a  run  with  capillarity.  The  same 
quadratic  saturation  profile  behind  the  5*  contour  (now  the  oil-water  front)  was  initialized  as 
occured  in  the  case  with  capillarity.  Fig.  6a  displays  H  and  FWHM  as  a  function  of  time  for 
the  5*  =  0.5  contour  for  such  a  run  with  M'  =  1.  It  is  superimposed  upon  the  analogous 
plots  for  the  5  =  0.0,  0.2,  0.4,  0.6,  0.8  contours  of  a  run  with  capillarity  P^  =  0.005.  Fig. 
6b  shows  the  same  superposition  for  M'  =  1.397.  The  effect  of  capillarity  is  found  to  affect 
the  rate  of  growth  of  finger  height  much  more  strongly  than  the  FWHM. 

4.2.   Stadies  in  Heterogeneous  Rock 

The  heterogeneity  properties  of  reservoir  rock  vary  according  to  rock  type.  We  employ 
here  a  characterization  of  heterogeneity  in  terms  of  a  length  scale  and  a  variance  of  the  rock 
properties  over  this  length  scale.    The  rock  heterogeneity  is  given  by  a  random  variable 
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(field)  K{x,  y)  =  K^^.  x{x,  y)  which  has  a  given  variance  about  the  mean  K^^  over  a  length 
scale  Z,;,,,.  A  'heterogeneity  mesh'  Ajc  =  A>  =  L^„  is  superimposed  on  the  computational 
square.  Uniformly  distributed  numbers  are  determined  at  the  nodes  of  this  mesh.  Each  uni- 
formly distributed  number  is  then  converted  to  a  normally  distributed  nimiber.  These  are 
then  interpolated  into  each  mesh  square  to  obtain  a  random  field  at  any  computational  point. 
Fig.  7  presents  a  plot  of  the  heterogeneity  field  produced  for  L;,„  =  0.05.  (The  numbers 
displayed  in  the  plot  are  relative.)  It  has  the  subjective  disadvantage  of  'not  looking  very 
heterogeneous',  especially  for  larger  /?;,„  values  (Compare  Figures  7  and  9b),  largely  due  to 
the  underlying  regular  grid  structure  and  linear  interpolation  by  which  it  is  constructed.  This 
formulation  does  however  isolate  a  single  length  scale  for  the  heterogeneity.  In  contrast  to 
real  rock  though,  any  interactions  due  to  interplay  between  different  length  scales  of  hetero- 
geneity eu'e  absent. 

The  variation  of  capillary  pressure  with  the  now  spatially  dependent  rock  permeability  is 
taken  into  account  using  the  Leverett  'J'-function  analysis.  Consequently,  for  the  studies  in 
heterogeneous  rock  Pq{x,  y)  in  (18)  becomes 

(.Xix,y))^^ 

where  P^^  is  a  constant. 

A  series  of  calculations  were  made  for  a  displacement  flow  of  given  frontal  mobility 
ratio  and  capillary  pressure  over  media  of  varying  heterogeneity  length  scale  for  constant 
heterogeneity  variance.  The  saturation  profile  was  initialized  as  a  one  dimensional  flow 
"o(^'  y)  ~  "oCy)  ^'^  ^  quadratic  variation  in  the  y  direction.  Three  runs  were  made  for  with 
M'  =  1.397  ,  P^^.  =  0.005  ,  with  a  very  strong  variance  in  the  heterogeneity  (permeability 
varying  by  a  factor  of  400) . 


run  Lhe,  I  ^ 

1  0.02 

2  0.05 

3  0.20 
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Figurcs  8a  plots  the  saturation  versus  injected  pore  volume  at  four  equally  spaced  time  inter- 
vals for  nm  1.  Figures  8b  and  8c  present  the  corresponding  plots  for  nms  2  and  3.  In  Fig. 
8a  the  capillary  length  scale  is  approximately  twice  as  large  as  the  heterogeneity  length  scale. 
In  spite  of  the  strong  variation  in  the  heterogeneity,  the  flow  is  stable,  the  variation  of  a  statura- 
tion  contour  does  not  grow  with  time.  In  Fig.  8b  the  heterogeneity  and  capillary  length  scales 
are  ajjproximately  equal.  In  Fig.  8c,  the  flow  is  cleary  unstable;  pockets  of  oil  are  being 
trapped  in  low  mobility  regions  as  neighbouring  fingers  flow  around  the  region  and  entrap. 
In  Fig.  9a  we  present  the  results  for  a  run  with  a  much  smaller  range  of  variation  variation  in 
the  heterogeneity.  For  this  ma  R^  =  0.0077,  R^p  a  0.03,  R^„=0.3.  and  the  heterogeneity 
(Fig.  9b)  in  the  permeability  varies  by  a  maximum  of  4.3  .  The  flow  is  still  dearly  unstable. 

5.   Conclusions 

We  have  presented  numerical  calculations  that  are  in  good  qualitative  agreement  with 
laboratory  studies  indicating  fingering  phenomenon  for  immiscible  flow  in  reservoir  rock 
samples.  Our  calcualtions  show  that,  although  immiscible  displacement  is  expected  to  be 
much  more  stable  than  misdble,  it  is  none-the-less  significantly  unstable  with  respect  to 
viscous  fingering  over  a  large  range  of  its  parameter  space.  Calculations  in  homogeneous 
rock  show  that  mobility  ratio,  capillary  pressure,  and  perturbation  size  are  factors  that  con- 
trol the  instability  of  the  flow.  While  growth  of  small  pertuirbations  is  largely  driven  by  the 
frontal  mobility  ratio,  the  growth  of  large  perturbations  cannot  in  general  be  characterized  by 
a  single  mobility  ratio,  though  the  frontal  mobility  ratio  tends  to  be  a  better  indicator  them 
the  far  field  mobility  ratio.  Both  imstable  growth  and  stable  decay  of  finger  height  are  found 
to  be  well  characterized  by  constant  rates. 

The  introduction  of  heterogeneity  in  the  medium  leads  to  the  evolution  of  finger  pat- 
terns with  the  proviso  that  fingers  induced  by  heterogeneity  on  length  scales  small  compared 
to  the  capillary  length  are  strongly  damped  by  capillarity;  the  medium  appears  effectivley 
homogeneous  over  these  small  length  scales.  Where  the  heterogeneity  length  scale  exceeds 
the  capillary  length  scale,  the  displacement  flow  'sees'  significant  mobility  channels  which  will 
act  as  nucleation  points  for  pertxirbations  in  the  displacement  flow.  The  subsequent  growth 
rate  of  these  perturbations  must  then  depend  upon  the  size  of  the  perturbation  achieved  as  a 


-15- 

resuJt  of  the  channel  and  the  mobility  ratio  of  the  fluids  involved. 

Comparison  with  nms  with  no  capillarity  show  the  expected  results  that  capillarity  tends 
to  broaden  existing  fingers  and  enhance  the  stability  of  the  flow. 

Nomenclature 

I  =    phase,  w(ater)  or  o(il) 

k^  =    relative  permeability,  phase  i 

k^  =   maximum  oil  relative  permeability  (S  =  0) 

k,^  =   maximum  water  relative  permeability  (5  =  1) 

K  =    rock  permeability 

X^  =    water  mobility 

X^  =    oil  mobility 

X  =    total  fluid  mobility 

L  =    computational  region  length  scale 

L      =    capillary  length  scale 

M  =    'far  field*  mobility  ratio 

L^„  =   heterogeneity  length  scale 

M'  =   frontal  mobility  ratio  defined  for  Buckley-Leverett  flow 

ii  =    viscosity 

o  =    oil 

Pg  -    oil  pressure 

P^  —   water  pressure 

P  -    mobility  weighted  average  pressure 

P^  =   capillary  pressure 

S„  -    oil  saturation 

5^  =    water  saturation 

Sy,c  —   connate  water  saturation 

5onv  =    residual  oil  saturation 

S  =    normalized  water  saturation  (0  ^  5  ^  1) 
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V  =    total  fluid  velocity 

V^  =    oil  velocity 

V^  =   water  velocity 

w  =   water 

X  =  heterogeneity  random  variable 
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Flgare  Captions 

Fig.  1  A  sketch  of  a  finger  illustrating  the  measurement  of  finger  height  and  full  width  at 
half  maximum  (FWHM). 

Fig.  2  a)  Contours  of  saturation  at  four  equally  spaced  timesteps  for  run  3  of  §4.1.3  .  b) 
Contours  of  pressure  (upper  left  hand  comer),  saturation  (upper  right)  and  velo- 
city magnitude  (lower)  for  the  final  timestep  of  run  3  of  54.1.3  .  c)  Cross  sec- 
tional profiles  along  the  mid-line  x  =  0.5  of  the  pressure  (upper  graph)  and 
saturation  (lower  graph)  for  5  equally  spaced  timesteps. 

Fig.  3  a)  Finger  height  versus  injected  pore  volume  for  three  mobility  nms  in  homogene- 
ous rock:  left  to  right  A/*  =  0.722,  1.0,  1.397.  For  each  run  the  contours 
5  =  0.8,  0.6,  0.4,  0.2,  0.0  are  represented  (respectively,  from  the  bottom  to  top  at 
r  =  0).  Plot  lines  that  terminate  early  in  time  indicate  the  contour  line  reaches 
y  =  L  shortly  thereafter.  The  vertical  and  horizontal  scales  are  the  same  for  all 
three  rtms.  b)  The  ratio  H  /  FWHM  versus  injected  pore  volume.  All  other  com- 
ments in  a)  apply. 

Fig.  4  a)  Finger  height  versus  injected  pore  volume  for  four  choices  of  initial  perturba- 
tion size  corresponding  (left  to  right)  respectively  to  the  runs  1  through  4 
described  in  §4.1.2  .  For  each  run  the  contours  S  =  0.8,  0.6,  0.4,  0.2,  0.0  are 
represented  (respectively,  from  the  bottom  to  top  at  f  =  0).  Plot  lines  that  ter- 
minate early  in  time  indicate  the  contour  line  reaches  y  =  L  shortly  thereafter. 
The  vertical  and  horizontal  scales  are  the  same  for  all  three  nms.  b)  The  ratio 
H  /  FWHM  versus  injected  pore  volume.   All  other  comments  in  a)  apply. 

Fig.  5  a)  Finger  height  versus  injected  pore  volume  for  five  choices  of  capillary  pressure 
strength  left  to  right  Pq  =  0.005,  0.01,  0.02,  0.05,  0.1  atm.  For  each  run  the  con- 
tours S  =  0.8,  0.6,  0.4,  0.2,  0.0  are  represented  (respectively,  from  the  bottom  to 
top  at  r  =  0).  Plot  lines  that  terminate  early  in  time  indicate  the  contour  line 
reaches  y  =  L  shortly  thereafter.  The  vertical  and  horizontal  scales  are  the  same 
for  all  three  nms.  b)  The  ratio  H  /  FWHM  versus  injected  pore  volume.  All 
other  comments  in  a)  apply. 
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Fig.  6  a)  Comparison  of  height  and  FWHM  versus  injected  pore  volume  for  the  5.  con- 
tour in  a  run  with  no  capillarity  with  the  analyzed  contours  for  a  corresponding 
run  with  Pq  =  0.005.  The  frontal  mobility  ratio  is  1.0  .  b)  As  in  a)  with  a  frontal 
mobility  ratio  of  1.397  . 

Fig.  7  The  heterogeneity  field  (contours  of  constant  heterogeneity)  with  a  length  scale 
Z,;,,,  /  Z,  =  0.05  ,  variance  =  1.0  . 

Fig.  8  a)  The  saturation  5(x,  y)  at  four  equally  spaced  timesteps  for  run  1  of  §  4.2  .  b) 
The  saturation  S(x,  y)  at  four  equally  spaced  timesteps  for  run  2  of  5  4.2  .  c)  The 
saturation  5(x,  y)  at  four  equally  spaced  timesteps  for  run  3  of  §  4.2  .  A  small 
amount  of  numerical  instability  has  produced  small  oscillation  behaviour  in  the 
trapped  oil  regions. 

Fig.  9  a)  The  saturation  S{x,  y)  at  four  equally  spaced  timesteps  for  the  final  nm  dis- 
cussed in  §  4.2  .    b)  The  heterogeneity  function  x{x,  y)  for  this  nm. 
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